home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / math / alged34.zip / ALGEDSRC.ZIP / ALGPOLY.C < prev    next >
C/C++ Source or Header  |  1996-06-06  |  16KB  |  540 lines
  1. /*--------------------------------------------------------------------
  2.    Alged:  Algebra Editor    henckel@vnet.ibm.com
  3.  
  4.    Copyright (c) 1994 John Henckel
  5.    Permission to use, copy, modify, distribute and sell this software
  6.    and its documentation for any purpose is hereby granted without fee,
  7.    provided that the above copyright notice appear in all copies.
  8. */
  9. #include "alged.h"
  10. /*-----------------------------------------------------------------
  11.    make-tree convert a table of coefficients to a node tree.
  12.    Note
  13.    the coefficients are REUSED (NOT copied), so don't free them.
  14. */
  15. node *maketree(node **coef, int sz, node *base) {
  16.   int i;
  17.   node *p,*tmp;
  18.  
  19.   p = coef[0];
  20.   for (i=1; i<sz; ++i) {
  21.     if (coef[i]->kind==NUM && !coef[i]->value) {
  22.       freenode(coef[i]);             /* throw away zeros */
  23.     }
  24.     else {
  25.       tmp = p;
  26.       p = newoper(ADD);
  27.       p->lf = tmp;
  28.       p->rt = tmp = newoper(MUL);
  29.       tmp->lf = coef[i];             /* add coef[i]*base^i */
  30.       if (i>1) {
  31.         tmp->rt = newoper(EXP);
  32.         tmp->rt->lf = deepcopy(base);
  33.         tmp->rt->rt = newnum(i);
  34.       }
  35.       else tmp->rt = deepcopy(base);
  36.     }
  37.   }
  38.   return p;
  39. }
  40.  
  41. /*-----------------------------------------------------------------
  42.    looks inside 'a' for any expression b or b^N.
  43.    if found, it changes it to a 1 (ONE) and returns the exponent.
  44.    Note, if N is not an integer (it's a fraction or expression)
  45.    then findbase ignores it.
  46.    When findbase returns 0, then  you can assume 'a' is unchanged.
  47. */
  48. int findbase(node *a, node *b) {
  49.   int i,x=1;
  50.  
  51.   if (equal(a,b) ||                         /* any expression */
  52.       a->kind==EXP && a->rt->kind==NUM &&       /* EXPONENTS */
  53.       (x = a->rt->value, whole(x)) &&
  54.       x >=0 && equal(a->lf,b)) {
  55.     movenode(a,newnum(1));
  56.     return x;
  57.   }
  58.   if (a->kind==MUL) {                     /* recurse on MULTIPLY */
  59.     x = findbase(a->lf,b);
  60.     if (x) return x;
  61.     x = findbase(a->rt,b);
  62.     return x;
  63.   }
  64.   return 0;
  65. }
  66.  
  67. /*-----------------------------------------------------------------
  68.    add entry to coef table          grow and init table
  69. */
  70. #define addcoef(k,a) do {                 \
  71.   node *tmp;                              \
  72.   if (sz<=i) {                            \
  73.     coef = realloc((void*)coef,(i+1)*sizeof*coef);   \
  74.     for ( ; sz<=i; ++sz)                  \
  75.       coef[sz] = newnum(0);               \
  76.   }                                       \
  77.   tmp = coef[i];                          \
  78.   coef[i] = newoper(k);                   \
  79.   coef[i]->lf = tmp;                      \
  80.   coef[i]->rt = a;                        \
  81. } while(0)
  82.  
  83. /*-----------------------------------------------------------------
  84.    MAKETABLE convert a tree to a table of coefficients which is indexed
  85.    by the power of a given base.  For example, a*x^2 + 3*x - b would be
  86.       0 -> 0 - b
  87.       1 -> 0 + 3
  88.       2 -> 0 + a
  89.    the size returned is the index of the last non-zero coefficient + 1
  90.    Note:
  91.    the tree, p, is expected to have correct association.
  92.    the coefficients always have a leading zero.
  93.    the tree is not altered or referred to by the table.
  94.    on return, you are guaranteed that result[i]->kind is one of NUM,
  95.       ADD, SUB and if it is NUM then it is zero.
  96. */
  97. node **maketable(node *base, node *p, int *size) {
  98.   node *a;
  99.   int i,sz;
  100.   node **coef;
  101.  
  102.   coef = malloc(sizeof*coef);        /* initial size = 1 */
  103.   checknull(coef);
  104.   *coef = newnum(0);
  105.   sz = 1;
  106.   a = deepcopy(p);                /* make a working copy */
  107.   while (1) {
  108.     i = findbase(a,base);
  109.     if (!i && aop(a->kind)) {
  110.       i = findbase(a->rt,base);
  111.       addcoef(a->kind,a->rt);
  112.     }
  113.     else {
  114.       addcoef(ADD,a);
  115.       break;
  116.     }
  117.     a = a->lf;
  118.   }
  119.   *size = sz;
  120.   return coef;
  121. }
  122.  
  123. /*-----------------------------------------------------------------
  124.    polynomial long division
  125.  
  126.    base is the expression on which the division is done (e.g. x)
  127.    nm and dn are numerator and denominator
  128.    returns:
  129.    a new tree if success, else returns NULL
  130. */
  131. node *polydiv(node *base, node *nm, node *dn) {
  132.   node **num,**den,**quo,**rem;
  133.   int szn,szd,szq,szr;
  134.   int i,j;
  135.   node *tmp,*p;
  136.  
  137.   num = maketable(base,nm,&szn);
  138.   den = maketable(base,dn,&szd);
  139.  
  140.   /*  If degree of numerator is less than denominator, or
  141.       if the denominator does not contain base, then return failure */
  142.  
  143.   if (szn < szd || szd < 2) {
  144.     for (i=0; i<szn; ++i) freetree(num[i]);
  145.     for (i=0; i<szd; ++i) freetree(den[i]);
  146.     free(num); free(den);
  147.     return NULL;
  148.   }
  149.   szq = szn-szd+1;
  150.   quo = malloc(szq*sizeof*quo);
  151.   checknull(quo);
  152.   szr = szd-1;
  153.   rem = malloc(szr*sizeof*rem);
  154.   checknull(rem);
  155.  
  156.   /*  Divide every coefficient by leading coeff in denominator */
  157.  
  158.   p = den[szd-1];
  159.   if (!(p->kind==ADD && p->rt->kind==NUM && p->rt->value==1.0)) {
  160.     for (i=0; i<szn; ++i) {
  161.       if (num[i]->kind != NUM) {      twirl();
  162.         tmp = num[i];
  163.         num[i] = newoper(DIV);
  164.         num[i]->lf = tmp;
  165.         num[i]->rt = deepcopy(p);
  166.       }
  167.     }
  168.     for (i=0; i<szd-1; ++i) {
  169.       if (den[i]->kind != NUM) {      twirl();
  170.         tmp = den[i];
  171.         den[i] = newoper(DIV);
  172.         den[i]->lf = tmp;
  173.         den[i]->rt = deepcopy(p);
  174.       }
  175.     }
  176.   }
  177.  
  178.   /*  Construct the coefficients of the quotient */
  179.  
  180.   for (i=1; i<=szq; ++i) {
  181.     quo[szq-i] = deepcopy(num[szn-i]);
  182.     for (j=1; j<i; ++j) if (i-j < szd) {      twirl();
  183.       tmp = quo[szq-i];
  184.       quo[szq-i] = newoper(SUB);
  185.       quo[szq-i]->lf = tmp;
  186.       quo[szq-i]->rt = tmp = newoper(MUL);
  187.       tmp->lf = deepcopy(den[szd-i+j-1]);
  188.       tmp->rt = deepcopy(quo[szq-j]);
  189.     }
  190.   }
  191.  
  192.   /*  Construct the coefficients of the remainder */
  193.  
  194.   for (i=0; i<szr; ++i) {
  195.     rem[i] = deepcopy(num[i]);
  196.     for (j=0; j<=i; ++j) if (i-j < szq) {      twirl();
  197.       tmp = rem[i];
  198.       rem[i] = newoper(SUB);
  199.       rem[i]->lf = tmp;
  200.       rem[i]->rt = tmp = newoper(MUL);
  201.       tmp->lf = deepcopy(den[j]);
  202.       tmp->rt = deepcopy(quo[i-j]);
  203.     }
  204.   }
  205.  
  206.   /*  Convert the quo and rem tables into a node tree */
  207.  
  208.   p = newoper(ADD);
  209.   p->rt = newoper(DIV);
  210.   p->rt->lf = maketree(rem,szr,base);
  211.   p->rt->rt = deepcopy(dn);
  212.   p->lf = maketree(quo,szq,base);
  213.  
  214.   /*  Calculate remainder trivial stuff */
  215.  
  216.   while (calcnode(p->rt,1));
  217.  
  218.   /*  Clean up and return */
  219.  
  220.   for (i=0; i<szn; ++i) freetree(num[i]);
  221.   for (i=0; i<szd; ++i) freetree(den[i]);
  222.   free(num);
  223.   free(quo);
  224.   free(rem);
  225.   free(den);
  226.   return p;
  227. }
  228.  
  229. /*-----------------------------------------------------------------
  230.    polycoef - collect the coefficients of a polynomial
  231. */
  232. void polycoef(node *base, node *p) {
  233.   node **coef;
  234.   int sz;
  235.   int i;
  236.  
  237.   if (!p || !base || p->kind==EQU || base->kind==EQU) return;
  238.  
  239.   coef = maketable(base,p,&sz);
  240.   if (sz < 2) {
  241.     if (sz) freetree(coef[0]);
  242.     free(coef);
  243.     return;
  244.   }
  245.   for (i=0; i<sz; ++i)              /* calculate the coefficients */
  246.     while (calcnode(coef[i],1));
  247.  
  248.   movenode(p,maketree(coef,sz,base));     /* replace p */
  249.   free(coef);
  250. }
  251.  
  252. /*-----------------------------------------------------------------
  253.    use quadratic equation to factor a degree-2 polynomial
  254.                                      2                       2
  255.      2                     b + sqrt(b - 4ac)       b - sqrt(b - 4ac)
  256.    ax + bx + c  ==>  a(x + -----------------) (x + -----------------)
  257.                                  2a                       2a
  258. */
  259. void quadratic(node *base,node *p) {
  260.   node **coef,*a;
  261.   int sz;
  262.   int i;
  263.  
  264.   coef = maketable(base,p,&sz);
  265.  
  266.   if (sz != 3) {
  267.     for (i=0; i<sz; ++i) freetree(coef[i]);
  268.     free(coef);
  269.     return;
  270.   }
  271.  
  272.   /* Construct the two binomials p = (ax + u)(x + v)  */
  273.  
  274.   a = newoper(ADD);
  275.   a->lf = deepcopy(coef[1]);                     /* b */
  276.   a->rt = newoper(EXP);
  277.   a->rt->lf = newoper(SUB);
  278.   a->rt->lf->lf = newoper(EXP);                  /* b^2 */
  279.   a->rt->lf->lf->lf = deepcopy(coef[1]);
  280.   a->rt->lf->lf->rt = newnum(2);
  281.   a->rt->lf->rt = newoper(MUL);                  /* 4ac */
  282.   a->rt->lf->rt->lf = newoper(MUL);
  283.   a->rt->lf->rt->lf->lf = newnum(4);
  284.   a->rt->lf->rt->lf->rt = deepcopy(coef[2]);
  285.   a->rt->lf->rt->rt = deepcopy(coef[0]);
  286.   a->rt->rt = newnum(0.5);
  287.  
  288.   movenode(p,newoper(MUL));
  289.   p->lf = newoper(ADD);
  290.   p->lf->lf = newoper(MUL);                      /* ax */
  291.   p->lf->lf->lf = deepcopy(coef[2]);
  292.   p->lf->lf->rt = deepcopy(base);
  293.   p->lf->rt = newoper(DIV);
  294.   p->lf->rt->lf = a;
  295.   p->lf->rt->rt = newnum(2);
  296.  
  297.   p->rt = newoper(ADD);
  298.   p->rt->lf = deepcopy(base);
  299.   p->rt->rt = newoper(DIV);
  300.   p->rt->rt->lf = deepcopy(a);
  301.   p->rt->rt->lf->kind = SUB;
  302.   p->rt->rt->rt = newoper(MUL);                  /* 2a */
  303.   p->rt->rt->rt->lf = newnum(2);
  304.   p->rt->rt->rt->rt = deepcopy(coef[2]);
  305.  
  306.   while(calcnode(p->lf->lf->lf,1));
  307.   while(calcnode(p->lf->rt,1));
  308.   while(calcnode(p->rt->rt,1));
  309.  
  310.   for (i=0; i<sz; ++i) freetree(coef[i]);
  311.   free(coef);
  312. }
  313.  
  314. /*-----------------------------------------------------------------
  315.    use cubic equation to factor a degree-3 polynomial
  316. */
  317. void cubic(node *base,node *pol) {
  318.   node **coef;
  319.   int sz;
  320.   int i,k=1;
  321.   node *a3,*a,*b,*c,*A,*B,*Q,*p,*q,*y1,*y2,*y3,*r;
  322.  
  323.   coef = maketable(base,pol,&sz);
  324.  
  325.   if (sz != 4) {
  326.     for (i=0; i<sz; ++i) freetree(coef[i]);
  327.     free(coef);
  328.     return;
  329.   }
  330.   a = cons(deepcopy(coef[2]),DIV,deepcopy(coef[3]));
  331.   b = cons(deepcopy(coef[1]),DIV,deepcopy(coef[3]));
  332.   c = cons(deepcopy(coef[0]),DIV,deepcopy(coef[3]));
  333.  
  334.   while(calcnode(a,k));
  335.   while(calcnode(b,k));
  336.   while(calcnode(c,k));
  337.  
  338.   a3= cons(deepcopy(a),DIV,newnum(3));
  339.   p = cons(deepcopy(b),SUB,cons(deepcopy(a),MUL,deepcopy(a3)));
  340.   q = cons(newnum(2),MUL,cons(deepcopy(a3),EXP,newnum(3)));
  341.   q = cons(q,SUB,cons(deepcopy(a3),MUL,deepcopy(b)));
  342.   q = cons(q,ADD,deepcopy(c));
  343.  
  344.   while(calcnode(p,k));
  345.   while(calcnode(q,k));
  346.  
  347.   Q = cons(cons(deepcopy(p),DIV,newnum(3)),EXP,newnum(3));
  348.   Q = cons(Q,ADD,cons(cons(deepcopy(q),DIV,newnum(2)),EXP,newnum(2)));
  349.  
  350.   while(calcnode(Q,k));
  351.  
  352.   if (Q->kind!=NUM || Q->value >= 0) {        /* not "irreducible" case */
  353.     A = cons(deepcopy(q),DIV,newnum(-2));
  354.     A = cons(A,ADD,cons(deepcopy(Q),EXP,newnum(0.5)));
  355.     A = cons(A,EXP,newnum(1.0/3.0));
  356.  
  357.     B = deepcopy(A);
  358.     B->lf->kind = SUB;
  359.  
  360.     while(calcnode(A,k));
  361.     while(calcnode(B,k));
  362.  
  363.     y1 = cons(deepcopy(A),ADD,deepcopy(B));
  364.     y2 = cons(deepcopy(y1),DIV,newnum(2));
  365.     y2->lf->kind = SUB;
  366.     y2 = cons(y2,MUL,cons(newvar("i"),MUL,cons(newnum(3),EXP,newnum(0.5))));
  367.     y2 = cons(cons(deepcopy(y1),DIV,newnum(-2)),ADD,y2);
  368.  
  369.     y3 = deepcopy(y2);
  370.     y3->kind = SUB;
  371.  
  372.     while(calcnode(y1,k));
  373.     while(calcnode(y2,k));
  374.     while(calcnode(y3,k));
  375.  
  376.     r = deepcopy(coef[3]);
  377.     r = cons(r,MUL,cons(deepcopy(base),SUB,cons(y1,SUB,deepcopy(a3))));
  378.     r = cons(r,MUL,cons(deepcopy(base),SUB,cons(y2,SUB,deepcopy(a3))));
  379.     r = cons(r,MUL,cons(deepcopy(base),SUB,cons(y3,SUB,deepcopy(a3))));
  380.  
  381.     movenode(pol,r);
  382.  
  383.     freetree(A);
  384.     freetree(B);
  385.   }
  386.  
  387.   for (i=0; i<sz; ++i) freetree(coef[i]);
  388.   freetree(a);
  389.   freetree(b);
  390.   freetree(c);
  391.   freetree(a3);
  392.   freetree(Q);
  393.   freetree(p);
  394.   freetree(q);
  395.   free(coef);
  396. }
  397.  
  398. /*-----------------------------------------------------------------
  399.    next factor, returns the ith factor of an expression.
  400.    it is a new nodetree.
  401. */
  402. node *nextfact(node *p,int i) {
  403.   int j;
  404.   double k;
  405.   if (i<1) return newnum(1);          /* everything has factor 1 */
  406.   if (p->kind==NUM && (k = fabs(p->value), whole(k))) {
  407.     if (k==1) return NULL;            /* don't return 1 twice */
  408.     for (j=2; j<=maxrat && j*2<=k; ++j)
  409.       if (fmod(k,j)==0 && !--i)
  410.         return newnum(j);
  411.   }
  412.   else if (p->kind==MUL) {
  413.     while (i>2 && p->kind==MUL) {
  414.       p=p->lf; i-=2;
  415.     }
  416.     if (i==2 && p->kind==MUL)
  417.       return deepcopy(p->rt);
  418.   }
  419.   if (i==1) return deepcopy(p);
  420.   return NULL;
  421. }
  422.  
  423. /*-----------------------------------------------------------------
  424.    check root
  425.    returns the new quotient on success, else null;
  426. */
  427. node *checkroot(node *base, node *p, node *den) {
  428.   node *ans;
  429.   twirl();
  430.   ans = polydiv(base,p,den);
  431.   if (!ans) return NULL;
  432.   if (ans->kind==ADD) {                /* always true!! */
  433.     while (distribute(ans->rt));
  434.     simplify(ans->rt);
  435.     while (calcnode(ans->rt,0));
  436.     simplify(ans->rt);
  437.     if (ans->rt->kind==NUM && ans->rt->value==0)    /* success!! */
  438.       return ans;
  439.   }
  440.   freetree(ans);
  441.   return NULL;
  442. }
  443.  
  444. /*-----------------------------------------------------------------
  445.    factor a polynomial using rational roots
  446.    Note: this doesn't work unless all the coefficients are
  447.    integers, or at least, not some form of quotient like (a*x/b).
  448. */
  449. void factrpoly(node *base,node *p) {
  450.   node **coef,*a,*b,*q,*dn,*nm,*r;
  451.   int sz;
  452.   int i,j,n;
  453.  
  454.   coef = maketable(base,p,&sz);
  455.   if (sz < 3) {
  456.     for (i=0; i<sz; ++i) freetree(coef[i]);
  457.     free(coef);
  458.     return;
  459.   }
  460.  
  461.   /* For every possible rational root, check it for divisibility  */
  462.  
  463. //  debug(coef[0]); debug(coef[sz-1]);
  464.   if (coef[0]->kind!=NUM) {
  465.     if (coef[0]->rt->kind==DIV &&  // coef is (0 + a/b) rational number
  466.         (a = coef[0]->rt->lf)->kind == NUM &&
  467.         (b = coef[0]->rt->rt)->kind == NUM &&
  468.         coef[sz-1]->kind!=NUM) {
  469.       coef[sz-1]->kind = MUL;
  470.       coef[sz-1]->lf = b;       // move b to the leading coef
  471.       if (coef[0]->kind==SUB)
  472.         a->value = -a->value;
  473.       freenode(coef[0]->lf);
  474.       freenode(coef[0]->rt);
  475.       freenode(coef[0]);
  476.       coef[0] = a;
  477.     }
  478.     else {
  479.       while (calcnode(coef[0],0));
  480.       if (coef[0]->kind!=NUM) {
  481. //        primefact(coef[0]);     this may cause trouble with exponents
  482. //        while (distribute_c(coef[0]));   // y^(2*2) = y^2 * y^2
  483.         while (exexpand(coef[0]));
  484.         while (nosubdiv(coef[0]));
  485.         while (distribute_c(coef[0]));   // (xy)^-1 = x^-1 * y^-1
  486.         while (fixassoc(coef[0]));
  487.       }
  488.     }
  489.   }
  490.   while (calcnode(coef[sz-1],0));
  491. //  debug(coef[0]); debug(coef[sz-1]);
  492.   nm = deepcopy(p);        /* initialize numerator */
  493.   n = 2;                   /* n is factor counter */
  494.   r = NULL;                /* intialize the list of factors */
  495.   for (i=0; n<sz; ++i) {
  496.     b = nextfact(coef[sz-1],i);     /* get next factor */
  497.     if (!b) break;
  498.     for (j=0; n<sz;) {
  499.       a = nextfact(coef[0],j);         /* next factor */
  500.       if (!a) break;
  501.       /*-----------------------------------------------------------------
  502.          check roots a/b, -a/b, ia/b, -ia/b
  503.       */
  504.       dn = cons(deepcopy(base),SUB,cons(a,DIV,deepcopy(b)));
  505.       simplify(dn->rt);                /* twice for good luck */
  506. #define DO_CHECK          \
  507.       simplify(dn->rt);    \
  508.       q = checkroot(base,nm,dn); \
  509.       if (q) {             \
  510.         r = cons(r,MUL,dn);\
  511.         freetree(nm);      \
  512.         nm = q;            \
  513.         while(distribute(nm)); \
  514.         simplify(nm);      \
  515.         simplify(nm);      \
  516.         ++n;               \
  517.         continue;          \
  518.       }
  519.       DO_CHECK
  520.       dn->rt = cons(dn->rt,MUL,newnum(-1));
  521.       DO_CHECK
  522.       dn->rt = cons(dn->rt,MUL,newvar("i"));
  523.       DO_CHECK
  524.       dn->rt = cons(dn->rt,MUL,newnum(-1));
  525.       DO_CHECK
  526.       freetree(dn);
  527.       ++j;
  528.     }
  529.     freetree(b);
  530.   }
  531.   if (n > 2) {
  532.     r = cons(r,MUL,nm);
  533.     movenode(p,r);           /* replace p with r */
  534.   }
  535.   for (i=0; i<sz; ++i) freetree(coef[i]);
  536.   free(coef);
  537.   return;
  538. }
  539.  
  540.